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Abstract 

To describe the strongly nonlinear dynamics of waves propagating in the final 
stages of shoaling and in the surf and swash zones, fully nonlinear models 
are required. The ability of the Serre or Green Naghdi (S-GN) equations to 
reproduce this nonlinear processes is reviewed. Two high-order methods for 
solving S-GN equations, based on Finite Volume approaches, are presented. 
The first one is based on a quasi-conservative form of the S-GN equations, 
and the second on a hybrid Finite Volume/Finite Difference method. We 
show the ability of these two approaches to accurately simulate nonlinear 
shoaling, breaking and runup processes. 
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1. Introduction 



Wave propagation in shallow water, and associated processes such as 
wave-breaking and run-up, play an important role in the nearshore dy- 
namics. The classical and successful method of describing slowly evolving 
wave-induced circulation in the nearshore is based on the phase- averaged 
approach, in which the depth-integrated mass and momentum equations are 
time-averaged over a wave period (see Phillips 49|). However, very important 
unsteady processes, such as wave run-up in the swash zone, coastal flooding 
during storm, tsunami and tidal bore propagation, require phase-resolving 
models. For coastal applications, these models are based on nonlinear shallow 
water equations (NSWE) and Boussinesq-type equations (see Brocchini and 
Dodd IJ])- NSWE give a good description of the nonlinear non-dispersive 
transformation of broken-waves, represented as shocks, in the inner surf and 
swash zones. However, due to the absence of frequency dispersion, the NSWE 
can not be applied to wave propagation before breaking. On the other hand, 
Boussinesq equations incorporate frequency dispersion and can be applied to 
wave shoaling, but contrary to NSWE they do not implicitly take into ac- 
count wave breaking. Since the 1990's, significant efforts have been devoted 
to extend the validity range of Boussinesq equations, by developing wave 
breaking parametrizations and by improving dispersive properties of these 
equations. Most of Boussinesq models used for nearshore applications are 
based on classical assumptions of weak nonlinearity, e = a/ho (a of the 
order of free surface amplitude, the characteristic water depth) and bal- 
ance between dispersion and nonlinearity: e = 0(/i) ^ 1, where /i = {ho/ L)"^ 
( L the characteristic horizontal scale). However, these assumptions may 
severely restrict applicability to real nearshore applications. Indeed, in the 
final stages of shoaling or in the surf zone, the wave dynamics is strongly 
nonlinear: e = 0(1). For instance, e is close to 0.4 just before breaking and 
can be larger than 1 in the swash zone. 

In 1953, a breakthrough treating nonlinearity was made by Serre (see 
Barthelemy (sj for a review). He derived ID fully nonlinear (e = 0(1)) 
weakly dispersive equations for horizontal bottom. Green and Naghdi [27| 
derived 2D fully nonlinear weakly dispersive equations for uneven bottom 
which represent a two-dimensional extension of Serre equations. Except for 
being formulated in terms of the velocity vector at an arbitrary z level, the 



equations of Wei et al. [60|] are basically equivalent to the 2D Serre or Green- 



Naghdi equations. It is now recognized that the Serre or Green Naghdi 
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(S-GN) equations represent the relevant system to model highly nonlinear 
weakly dispersive waves propagating in shallow water (see Lannes and Bon- 
neton j36|). However, much remains to be done for a proper representation 
of wave breaking, and for an accurate modeling of moving shorelines over 
strongly varying topographies. 

Due to facility of implementation, most of the Boussinesq-type models 
use finite difference schemes to discretize the equations (e.g. Abbott et al. 
[H , Wei et al. 60|). Finite volume methods, which are derived on the basis 
of the integral form of the conservation law, have many advantages. They 
are conservative and easily formulated to allow for unstructured meshes. Al- 
though the finite volume method has been widely jLnd successfully used to 
solve the strictly hyperbolic NSWE 
Marche et al. 
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4|) 



e.g. Leveque (38|, Brocchini and Dodd 



its application to Boussinesq-type equations has 

Stansby |55|). This 



only been reported recently (Bradford and Sanders [11 



is related to the fact that finite volume methods essentially aim at a good 
representation of advection, while methods used for Boussinesq-type equa- 
tions must also deal with third order derivatives responsible for dispersive 
effects. To overcome this problem, hybrid approaches, coupling finite vol- 
ume and the finite difference methods have been recently proposed (Soares 



Frazao and Zech 54 



25|). The hyperbolic 



Bernetti et al. [6|, Erduran et al. 
terms are treated using shock- capturing methods while the dispersive terms 
are discretized using the finite difference formulation. Until now, in coastal 
applications finite volume and hybrid approaches have only been applied for 
weakly nonlinear forms of Boussinesq-type equations. However, in the final 
stages of shoaling and in the surf and swash zones, the effects of nonlinearity 
are too large to be treated as a small perturbation. In the context of the 
IDAO ocean research programme we have investigated, the last few years, ap- 
plications of finite volume and hybrid methods to the fully nonlinear weakly 
dispersive S-GN equations. In this paper, we present a synthesis of this work, 
emphasizing the ability of S-GN models to deal with wave transformation in 
the surf and swash zones. 



2. Theoretical background 

According to jij and 36|, the 2D S-GN equations can be written under 
the following non-dimensionalized form 

Ct + V-(/iv) = 

vt + £(v v)v + vc = -^lV , (1) 
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where Ci'^yt) is the surface elevation, /i(x, t) = 1 + e( — b the water depth, 
6(x) the variation of the bottom topography and v(x, t) = {u,v) the depth 
averaged velocity. V characterizes non-hydrostatic and dispersive effects and 
writes 



V = r[h, b]v, + e (^-i-V {h%{^r . V)(V ■ v) - (V ■ v)^)) + Q[h, 6](v)^ 



(2) 



where the linear operator T[h, b] is defined as 



r[h, b]W = -^V{h^V-W) + ^[V{h^Vb-W)-h'^VbV-W] +VbVb-W (3) 
and the purely topographical term Q[h, 6](v) is given by 

Q[/i,&](v) = i-[v(/i2(vV)26) -/i2((v V)(V-v)-(V-v)2)V6] 

+ ((v- V)26)V6. (4) 

The range of validity of this set of equations can be easily extended to 
wave propagation problems in deeper waters using the dispersion correction 



technique discussed in references [6l[, [4l[, [18| and [10|. The frequency dis- 



persion properties are improved by applying the operator I + fi{a — l)T[h, b] 
to the momentum equation ([1]) and neglecting 0(yU.^) terms, which makes the 
term, 

(a-l)/i^ = -(a-l)/ir[/i,6](vt + £(v- V)v + VC) , (5) 

appears on the right-hand side of the momentum equation ([T]). 
The coefficient a is an adjustable parameter that must be tuned in order 
to minimize the phase and group velocity errors in comparison with the 
linear Stokes theory. This yields the optimal value a = 1.159 0. Inspired by 



Nwogu |46|, it is also possible to choose another dependent variable (such 



as the velocity at a certain depth) rather than the mean velocity, as in [6 



ITj . This allows for more freedom to match the Stokes linear dispersion 
and/or the exact linear shoaling coefficient. The price to pay is that the 
mass conservation equation is not exact anymore, but of order 0(/i^). 

It is known that the S-GN equations are mathematically well-posed in the 
sense that they admit solutions over the relevant time scale for any initial 



"'^the parameter used in reference [20| is given by a' = (q — l)/3, with an optimal value 
a' = 0.053 
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data reasonably smooth (see Isl for the general case, and simpler proofs 
for one dimensional surfaces [39|, |30|). Moreover, the solution of the S-GN 
equations provides a good approximation of the solution of the full water 
waves equations (see [2| for the general case, and 39| for one dimensional 
surface waves over fiat bottoms); this means that the difference between 
both solutions remains of order O(yU^) as long as the wave does not exhibit 
any kind of singularity such as wave breaking. Near the breaking point, 
the relevance of the S-GN equations is a completely open problem since the 
approximations made to derive the non-hydrostatic and dispersive terms ([2]) 
may diverge. Comparison of numerical simulations with experimental data 
are therefore necessary to assess the validity of the S-GN equations near wave 
breaking. 

The nonlinear shallow water equations (NSWE), or Saint Venant equa- 
tions, are obtained when the dispersive term /iP is neglected in the S-GN 
equations. It is well known that these nonlinear hyperbolic equations result 
in discontinuous solutions (shocks), which can be considered as the math- 
ematical counterparts of breaking wave front. Based on this idea, Hibbert 



and Peregrine [29| numerically simulated the entire process of bore propa- 



gation and runup on a constant beach slope. Kobayashi et al. [33|] applied 
the same approach to simulate the propagation of periodic broken-waves in 
the surf zone and found good results in comparison with laboratory data. A 
detailed analysis of the ability of the ID NSWE shock-wave model to pre- 
dict cross-shore wave transformation and energy dissipation in the inner surf 
zone was presented by Bonneton j^. For 2D problems. Peregrine 47] showed 
that non-uniformities along the breaking wave front (i.e. along the shock) 
due to alongshore inhomogeneities in the incident wave field or in the local 
bathymetry, drive vertical vorticity. Biihler 15|] presented a general theoret- 
ical analysis of wave-driven currents and vortex dynamics due to dissipating 
waves. From computations and laboratory measurements, Brocchini et al. 



13l | and Kennedy et al [32| showed that breaking- wave-generated vortices are 



qualitatively well described by the NSWE shock-wave theory. Bonneton et 
al. emphasized the importance of alongshore inhomogeneities of breaking 
wave energy dissipation for wave-induced rip current circulation. Their anal- 
ysis was based on the derivation of an equation for the mean-current vorticity, 
where the main driving term is related to shock-wave energy dissipation. 

The description of shallow water wave dynamics in realistic situations, 
i.e. over uneven bathymetries from the shoaling zone up to the shoreline, 
requires the development of advanced numerical approaches to integrate fully 
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nonlinear Boussinesq-type equations. In this framework, the S-GN equations 
offers the quite exceptional property of admiting closed form solutions of 
solitary and cnoidal type, bringing the opportunity to assess the accuracy 
and efficiency of numerical methods. 

For horizontal bottoms, the ID S-GN equations have an exact solitary- 
wave solution given by, in dimensional variables. 



h{x,t) = ho + Hsech^{K{x — ct)), 



u{x, t) = c( 1 



ho 



h{x, t) 



K 



2hoVho + H'' 

Vaih + H), 



(6a) 
(6b) 

(6c) 
(6d) 



where ho denotes the mean water depth and H the wave height. This family 
of solutions is known as the Rayleigh solitary wave solution 50|. Guizien 



and Barthelemy [28|] experimentally checked that solitary waves generated 
according to Rayleigh's law display very little dispersive trailing waves com- 
pared to KdV ones for instance. 



El et al. 2J] and Carter and Cienfuegos 16|] have recently shown that 
the S-GN equations admit also the following family of periodic solutions. 



h{x, t) = ao + aidn^ (k(x — ct), k) , 



u{x, t) = cil 



ho 



h{x, t) 



K 



^3ai 



2^Jao{ao + ai)(ao + (1 - k'^)ai) 
^^gaoiao + ai)(ao + (1 - k'^)ai) 



hn 



(7a) 
(7b) 



(7c) 



(7d) 



where k G [0, 1], Qq > 0, ai > are real parameters. In equation ([7]) dn(-, k) 
is a Jacobi elliptic function with elliptic modulus k. This family of solutions 
constitutes an important extension of the classic KDV cnoidal theory to 
strongly nonlinear and weakly dispersive applications. It is useful to relate 
the parameters of this solution to physical variables in order to compute 
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cnoidal waves in terms of wave height, H, wave period, T, and mean water 
depth, Hq. The latter is achieved by solving the following system of equations 
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«o = 'io-ai-^^ (8b) 

uj = K (8c) 

4[aoK{k) + aiE{k)f ' 

where Gj = Itx jT is the angular frequency, K{}i) and E(}z) are the complete 
elliptic integrals of the first and second kinds respectively. 

As /c — 1~, the family of periodic solutions limits to the two-parameter 
family of solitary- wave solutions given by equations (|6^-d). 

3. Reformulations of S-GN equations for numerical implementa- 
tions 

The formulation of the S-GN equations in function of conventional un- 
knowns {h,\), system ([1]), is not suitable for finite volume methods. In this 
section, we present two other formulations which are convenient for these 
numerical methods. 

S-GN equations in a quasi- conservative form 

In the ID case, it is possible to show that continuity and momentum equa- 
tions can be recast in a weak quasi-conservative form by defining an auxiliary 
variable q which aggregates all time derivatives in the momentum equations 



of system ([T]) (Cienfuegos et al. 18|). This convenient mathematical form 
writes down as, 

ht + F, = 0, 

qt + G,, = a'fiS', (9) 

where the source term in the right hand side is related to the small dispersive 
correction term presented in section 121 

S' = -2e{b - l)bA {uu, + . (10) 
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The auxiliary variable reads, 



u (11) 



g = <; 1 + /i + hi + hih^x - hhA - + - 1)^^ 

and the functions F and G are defined as follows, 

F = ehu, 

G = £{qu + C-^u^) 

It is important to note that if no dispersion correction is considered (i.e. 
a' = 0) S-GN equations can be written in the form of conservations law even 
if bottom variations are allowed. 

In the framework of numerical modelling, the system (|9]) can be conve- 
niently integrated over control volumes. 

S-GN equations in terms of the {h,h\) variables 

An alternative approach proposed by [lo[ is to write the S-GN equations 
in terms of the conservative variables [h, h\), namely 

ht + eV ■ (/iv) = 

(/iv)t + eV ■ (/iv ® v) + WC = -{I + fiahr[h,b]y)-'\-hVC 

II a 

+£/i/iQi[/i,6](v)l +-WC (12) 

a 

with Qi[/2,,6](v) = 6](v) — T[/i,6]((v ■ V)v). It is worth pointing out 
that this formulation does not include any third-order derivative, allowing for 
easier and robust numerical computations, especially when the wave becomes 
steeper. Note that the S-GN equations with improved dispersion a la Nwogu 



can also be put under a similar form [17 



This formulation is well-suited for a splitting approach with a finite vol- 
ume method for the hyperbolic part of the equations (the left-hand side of 
system (fT2|) ) and finite difference method for the dispersive part (the right- 
hand side of system f fT2|) ). 
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4. High-order compact finite volume method 



In the ID case, the quasi-conservative form iQ can be numerically inte- 
grated to describe wave propagation in shallow waters. However, in order 
to extend the application of the model into the surf zone, additional terms 
have to be added to the mass and momentum conservation equations. These 
terms aim at modelling wave-breaking energy dissipation and bottom fric- 
tion. In dimensional variables, this extended system can be written in the 
following form. 



qt + 



yDhu 
h 



ph 



a D , 



(13) 



where p is the water density, and D^u represent breaking terms, is the 
bed shear stress. 

Breaking-induced energy dissipation mechanisms are introduced through 
diffusive-like terms, Dh and Dhu, applied locally on the wave front face where 
an explicit breaking criterion is required to switch them on. The mathemat- 
ical form for Dh and Dhu is chosen in order to ensure that the overall mass 
and momentum budget is preserved, acting only as to locally redistribute 
these quantities under the breaker 2l|. Breaking terms are thus written in 
the form, 

Dh = {vhhx) , 
Dhu = dx {vhu{hu)x) , 

where Vh and Vhu are diffusivity functions expressed as. 



Vh{X) = -Khexp i- 1 



X 

Uhu{X) = -Khuexp {- 1 



X 



1 + 



X 



X 



X 

T 



with Kh and Khu slowing varying scaling coefficients, X is a moving hori- 
zontal coordinate attached to the wave crest and is the extent over which 
breaking terms are active. This breaking model has been calibrated on the 
Ting and Kirby's 59|] regular wave experiment and optimal parameter val- 
ues suggested by Cienfuegos et al. |2l| are Kh = 2cd, Khu = 20cd and 
Ir/d = 0.82, with c = (gd)^-^ and d the local still water depth. 
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The numerical integration of the system fll3l) is performed using a high- 
order compact finite volume method. A detailed description of this model, 



SERR-ID, is given in references [18|, |20[. The equations are first integrated 



in space over discrete control volumes fli = {x E [x^.i , 



di 



X. , 1 PX. , 1 



/idx + , t) — , t) = / dx{jyhhx)dx, (14) 



X . I J X . I 



— J qdx + G{Xi^i,t) - G{Xi^i,t) = J \^S' + -dx{i^hu{hu).x) - — ] dx, 



X. I J X . \ 



(15) 



where integral of variables h and q over control volumes must be advanced 
in time. The average value of function h at time t = tn over control volume 
is noted as, 



1 r^+h 

= ^ / ^(^' ^n) dx. 



where Ax = x,- , i — x,_i is the length of the discrete control volumes. Hence, 
integrated over the whole domain, equations (IT^ and (IT^ can be expressed 

as. 



dh? 1 



(i^+i - - {i^khx}Ui + {^hhx}U-l 



dt Ax 
dn^ ^ 1 / \ 

where N is the total number of control volumes used to discretize the physical 
domain and is the discretized counterpart of the source term in the right 
hand side of equation f|T5|) . This term is approximated through centred finite 
differences. The values h and q at cell interfaces are reconstructed from cell- 
averaged values using the irnplicit 4th order compact interpolation technique 



described in references |34l . l35| . At each time step, the velocity component 
at control volume interfaces is computed numerically by inverting equation 
(ITT]) . Time stepping is performed through a 4th order Runge-Kutta method. 

Efficient absorbing-generating boundary conditions have been implemented 
in SERR-ID. They are based on the following characteristic form of the S-GN 
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equations, 



^ = -^^(^'^)-^^^-;^ f = "-^' 

with positive and negative Riemann variables defined respectively as = 
u + 2y^^ and R~ = u — 2y/gh. Vertical acceleration of fluid particles, 
which is of order 0(/i) and thus disregarded in the nonlinear shallow water 
equations (NSWE), is represented by function P in the first term of the right 
hand side of equations f|T6|) - f|T7|l . This term is responsible for the loss of 
hyperbolicity in the S-GN equations by introducing an horizontal dependence 
in the characteristic plane {x,t). 

It is worth noting that from a physical point of view we can expect that 
time scales associated to dispersive effects would be larger than the ones 
associated to nonlinearities from intermediate to shallow waters. We may 
therefore assume that over short distances/times, Riemann variables might 
be locally conserved along characteristics. This physical argument has been 
used to develop absorbing-generating boundary conditions for the numerical 



resolution of S-GN equations 20|, |45 



For the moving shoreline boundary condition, the simple extrapolation 



technique proposed by Lynett et al. j40| has been adapted in the finite 
volume resolution. 

SERR-ID has been extensively validated by comparisons with non-breakinj 



and breaking wave laboratory experiments [19|, |20|, |21|, |45 . 

In the present paper, the capabilities of the model are illustrated by 
comparing numerical computations with breaking random wave propagation 
experiments. We use measurements conducted at the 70-meter long wave 
tank of the Instituto Nacional de Hidraulica (Chile), prepared with a beach 
of very mild slope of 1/80 in order to produce large surf zone extensions (23l |. 
A random JONSWAP type wave field (h0=0.52m, fp=0.25Hz, Hmo=0.17m) 
was generated by a piston wave-maker and measurements of the free surface 
displacements were performed all over its length at high spatial resolution 
(0.2m to Im). This experiment allows us to test numerical models that are 
assumed to reproduce nonlinear shallow water wave propagation, breaking 
and run-up. The evolution of the wave energy power spectral density as 
the wave field propagates over the beach is presented in Figure [H for both 
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x=-22[m] 



x=-9[m] 





x=9[m] 




Figure 1: Random wave transformation and breaking over a mild slope beach. The origin 
of the x-coordinate is at the toe of the beach slope, positive towards the shore. Blue : 
experimental data. Red : computed results. 



experimental measurements and numerical results. SERR-ID is able to sim- 
ulate the complex nonlinear energy transfer occuring in the shoaling and surf 
zones. In particular, it reproduces the generation of higher frequency har- 
monics during shoaling (between x=-9m and x=9m), energy dissipation by 
breaking in the surf zone (between x=9m and x=32m). The numerical results 
also indicate that the model is able to reproduce the energy transfer from 
the Jonswap spectrum band (>0.1Hz) to the infragravity band (<0.1Hz). It 
is important to note that the numerical model was forced with the high-pass 
filtered wave signal measured two meter away from the wave paddle without 
energy content in the infragravity band. 

5. High-order hybrid finite volume / finite difference method 

The formulation of S-GN equations introduced in section [3] (equations 
f lT2|) ) is well-suited for a splitting approach separating the hyperbolic and the 
dispersive part of the equations. In this section, we first present an efficient 
high-order positive preserving well-balanced shock-capturing scheme for the 
hyperbolic step and then the splitting method for solving the whole S-GN 
system. 
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NSWE shock capturing solver 

We consider in this section the hyperbohc part of the S-GN equation, 
stated in the dimensionnahzed form: 

ht + V ■ (/iv) = (18) 
{hv)t + V ■ (/iv ® v) + ghVC = 0. 

This system can be also regarded as an hyperbohc system of conservation 
laws with a geometrical source term controlled by the topography variations. 
To simplify the algorithm presentation we only consider the one- dimensional 
problem: 

wt + /(w), = ^(w), (19) 

where w = * (/i, hu), /(w) = * {hu, hv? + Q.hgh?) and S" = * (0, —ghhx) the 
source term. To perform numerical approximations of the weak solutions 
of this system, we use a high order finite-volume approach in conservative 



variables, relying on Riemann problems for hyperbolic conservative laws [26 
This approach allows accurate computation of propagating bores, with re- 
duced spurious effects of numerical dissipation and dispersion. Using such 
accurate scheme, we are able to handle wave breaking (see Section [2]). Since 
we aim at computing the complex interactions between propagating waves 
and topography (including the preservation of motionless steady states), we 
also embed this approach into a well-balanced scheme. 

More precisely, based on discrete finite- volume cell averaging w" at time 
t"' = n5t, we use the limited 4*'^-order MUSCL reconstruction suggested in 
0]. Considering a cell Ci, this approach provides, for all t", high order 
accuracy interpolated quantities Wj ^ and Wj ,., respectively at the left and 
right boundary of each cell. To get a positive preserving and well-balanced 
scheme, additional faces reconstructions are introduced j^: 

h* = max{bi^r,bi+i^i), 

h*,. = max(0, hi^r + Kr - b*), 

h*^^ i = max(0, hi+i^i + bi+i^i - b*). 

These new left and right values for water height are used to compute auxiliary 
conservative faces values w*^ and ^i+n- 
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which are injected into a Riemann solver. Neglecting temporally the time 
discretization, we obtain the following semi-discrete finite-volume scheme for 



where f - , i and f,_i are the numerical fiux functions based both on a conser- 
vative fiux consistent with the homogeneous NSWE issued from a relaxation 
approach [3], and the hydrostatic reconstruction correction to the interface 
fiuxes j^. Sc,i is a centered discretization of the source term needed to achieve 
accuracy, consistency and well-balancing properties. 

The resulting finite- volume scheme provides high-order accuracy approx- 
imations of the weak solutions of system f lT^ while preserving the positivity 
of the water-height, thanks to the relaxation approach developed in (3]. This 
last property is paramount to ensure robust and accurate simulations of time- 



evolving shorelines [4^. In addition, the use of a well-balanced scheme leads 
to accurate computations of incident waves and topography complex inter- 
actions classically occurring in the surf zone j43|. It also allows far more 
accurate results in situations involving tiny oscillations near steady states. 
The resulting high-order positive preserving well-balanced shock-capturing 
scheme was implemented in the code SURF-WB 44, 0]. 

We can observe on Figure [2] a comparison between numerical results and 
analytical solutions for a one dimensional test case. This solution describes 
time oscillations of a forced fiow over a quadratic topography profile, involv- 



ing a moving shoreline [51j. The free surface is always planar during the 
oscillations. The computational domain is 4320 m long and the topography 
is given by 

b{x) = 1 - /lo (- 



where Hq = 10 m and a = 3000 m are respectively the mean water depth 
and a topography scaling parameter. The fiow motion is driven by the left 
boundary condition, where we impose a periodic motion: 



'^--Wh-A^ cos(2v/8^*) ) - - (21) 
where i? is a free parameter, set to 2m.s~^ for the simulation shown here. 
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Note that the shorehne location, be analytically derived: 




( 



B^/Sgho/a"^ cos(^/8gho/aH) ] +a 



The simulation has been performed with 200 cells and the time step is set 
to 0.03 s, with a first order scheme. Results are shown first as a compar- 
ison between numerical results and analytical solution for the free surface 
elevation, at several times. Note that the results and the solution are almost 
indistinguishable. Then we highlight the accuracy of the shoreline location 
prediction through a comparison between analytical and numerically pre- 
dicted shoreline location, during 3000 s. 



Figure 2: (a) Free surface at various time between a half evolution period (analytical 
solution in solid line) (b) Shoreline location between i = s and t = 3000 s 

S-GN splitting solver 

Equations f lT2|) are solved using the following splitting method : we de- 
compose the solution operator S{-) associated to f|T2|) at each time step by 
the second order splitting scheme 



where Si{t) is the solution operator associated to the nonlinear shallow water 
equations ( |T8|) . while 5*2 (t) is the solution operator associated to the disper- 
sive part of the equations, namely 





(22) 
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As described previously, Si(t) is computed using a finite- volume ap- 
proach, while a finite-difference approach is used to solve 5*2 (t) at each time 
step: the spatial derivatives are discretized using centered fourth-order for- 
mulae. Boundary conditions are imposed using the method presented in (lo| . 
As far as time discretization is concerned, we choose to solve Si{t) and 5*2 (t) 
using an explicit fourth-order Runge-Kutta scheme. 

To sum up, Si(t) and 5*2 (t) are solved, within our splitting approach, using 
a fourth-order scheme in space and time. However, the use of a second-order 
splitting method implies that the global scheme is of order two in time. The 
use of a fourth-order Runge-Kutta scheme for Si{t) and 5*2 (t) is however 



required to have a good semi-discrete dispersion relation (see |10|). 

In order to handle wave breaking, we switch from the S-GN equations to 
the NSWE, locally in time and space, by skipping the dispersive step 5*2 
when the wave is ready to break. In this way, we only solve the hyperbolic 
part of the equations for the wave fronts, and the breaking wave dissipation 
is represented by the shock energy dissipation (see also |8|). To determine 
where to suppress the dispersive step at each time step, we use the first 
half-time step 5*1 of the time-splitting as a predictor to assess the local en- 
ergy dissipation. This dissipation is close to zero in regular wave regions, 
and forms a peak when shocks are appearing. We can then easily locate the 
eventual breaking wave fronts at each time step, and skip the dispersive step 
only at the wave fronts. 



In order to test the efficiency of our numerical methods, the ability of 
the model to describe the propagation of a strongly non-linear cnoidal wave 
solution of the S-GN equations (cf. relations ([7])) is investigated. Periodic 
boundary conditions have been used in order to appreciate the propagation 
of the cnoidal wave at long time, and the computational domain length is 
equal to the cnoidal wave-length. 

The initial condition is a cnoidal wave with H = 0.6m, ho = Im and 
T = 4s (see Figure [3]). Numerical and theoretical solutions are compared at 
t = 15T = 60s. Relative amplitude and celerity errors are given in Table [T] 
for different spatial steps. The Courant number remains equal to 1 for the 
different cases. 

Figure [3] shows the results for the finest grid size considered {dx = 0.01m). 
The cnoidal wave at t = Os and t = 60s are both plotted in this figure, 
but cannot be distinguished since the errors on wave height and celerity are 
extremely small (cf. Table [T]). In particular, the celerity error cannot be 
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precisely quantified for dx = 0.01m since tlie spatial lag between numerical 
and analytical solution at t = 60s is much smaller than the grid size : the 
numerical solution converge to the exact one for very small dx. Moreover, 
relative celerity and amplitude errors remain small for larger spatial steps, 
demonstrating the very high accuracy of our numerical methods. 

1.5 
1.4 
1.3 

1.2 
1.1 



0.9 



'-8 -6 -4 -2 2 4 6 8 

X (m) 

Figure 3: Propagation of a strongly nonlinear cnoidal wave over a periodic domain. At 
t = Os (blue line), the cnoidal wave is such as H = 0.6m, ho ^ Im and T = 4s. The red 
line represents the wave 15 periods later {t = 60s). The doted line is the still water level. 
dx ~ 0.01 TO and the Courant number is equal to 1. 






Relative Amplitude Error (%) 


Relative Celerity Error (%) 


dx = 0.01 m 


3.1 10-5 


< 5.4 10-3 


dx = 0.10 m 


-0.24 


0.026 


dx = 0.15 m 


-1.75 


0.031 



Table 1: Amplitude and celerity errors relative to the analytical solution for the cnoidal 
wave previously described (cf. Figure [S]) at t = 60s. Computations arc performed for 3 
spatial steps and a Courant number equal to 1. 

In the next case, we assess the ability of our model to describe wave 
runup and breaking. It is based on laboratory experiments carried out by 



Synolakis 58[, for an incident solitary wave of relative amplitude ao/ho 



0.28, propagating and breaking over a planar beach with a slope of 1:19.85. 
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The still water level in the horizontal part of the beach was = 0.3m, and 
the simulations are performed using the grid size dx = 0.08m and dt = 0.02s. 
As friction effects are important when the water becomes very shallow, in 
particular for the run-up and run-down stage, a quadratic friction term is 
introduced for this simulation. 

The comparison between measured and computed waves is presented in 
Figure |H It shows a good agreement between model predictions and labora- 
tory data for the wave shoaling, breaking, run-up and run-down. In partic- 
ular, the model is able to accurately describe the formation and breaking of 
the back- wash bore without any additional treatment, which is a particularly 
demanding test. 



6. Conclusion 

As the waves approach the shore, the nonlinearity effects become intense, 
especially in the final stages of shoaling and the surf zone. To simulate 
such nonlinear processes in shallow water, fully nonlinear Boussinesq-type 
approaches are required. The Serre or Green Naghdi (S-GN) equations, with 
improved dispersion properties, represent the relevant system to model these 



highly nonlinear weakly dispersive waves (see Lannes and Bonneton |36|). 
Two high-order methods for solving S-GN equations, based on Finite Volume 
approaches, are presented in this article. 

The first one is based on a quasi-conservative form of the S-GN equa- 
tions. Wave-breaking energy dissipation is taken into account through a 
diffusive-type parametrization on both the mass conservation and momen- 



tum conservation [21[. This model, SERR-ID, has been extensively validated 
by multiple comparisons between numerical simulations and physical exper- 
iments including solitary waves shoaling, regular waves propagating over a 



submerged bar [20( or regular wave breaking over uniform beach slopes [21 
In the present paper, new results show the ability of the model to reproduce 
nonlinear energy transfer for random waves, in the shoaling and surf zones. 

We present also an alternative approach where the S-GN equations are 
reformulated in terms of the conservative variables {h, hv) (equations (fT2|) ). 
This formulation is well-suited for a splitting approach with a finite volume 
method for the hyperbolic part of the S-GN equations (NSWE) and a fi- 



nite difference method for the dispersive part [10|. Our hyperbolic method 



is based on a high-order well-balanced shock-capturing scheme (SURF-WB 
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t*=10 



t*=20 




Figure 4: Comparisons of model predictions ( — ) and experimental data (x) for a breaking 
solitary wave with non-dimensional initial incident amplitude ao/^o = 0.28, on a 1 : 19.85 
constant slope beach investigated by Synolakis (1987). t* = t{g/hQy^^. Adapted from 
0. 
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Figure 5: Tidal bore propagation in Garonne River; Mascaret project (Parisot and Bon- 
neton, 2009-2011). 



code, jH, 0]). In the present article we show that our hybrid model accu- 
rately describes strongly nonlinear S-GN cnoidal wave solutions. In order to 
handle wave breaking, we switch locally from the S-GN equations to the hy- 
perbolic NSWE, where the wave is ready to break. In this way, we only solve 
the hyperbolic part of the equations for the wave fronts, and the breaking 
wave dissipation is represented by shock energy dissipation. We show that 
this approach accurately predicts nonlinear shoaling, breaking and runup 
of solitary waves on a beach. The advantage of this approach, in compar- 
ison with classical breaking parametrizations, is that it is easily extended 
to 2DH broken-wave problems. This is crucial to predict wave-induced cir- 
culations and macro- vortices, which are strongly controlled by dissipation 
non- uniformities along broken- wave fronts IJ, l9| . 

Further work is required to evaluate the capability of our S-GN models to 
predict such 2DH flows. An another important open problem concerns the 
ability of our approaches to predict bore dynamics in a large range of Froude 
numbers, from "non-breaking undular bore" to "breaking bore". The transi- 
tion between those two types of bores, which is controlled by the competition 
between nonlinearity, dispersion and dissipation effects, is still poorly under- 
stood. Recent field experiments on the dynamics of tidal bores (see figure [5]) 
will give us the opportunity to validate our models. 
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